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ABSTRACT: We construct models of the inner part of a transonic adiabatic accretion 
disc assuming constant specific angular momentum taking the vertical structure fully into 
account. 

For comparison purposes, we construct the corresponding one dimensional viscous disc 
1^ -_ models derived under vertical averaging assumptions. The conditions under which a unique 

j^ ■ location for the critical/sonic point is obtained, given an appropriate set of exterior bound- 

ary conditions for these models, is also discussed. This is not unique if the standard 'a' 
prescription with viscous stress proportional to the angular velocity gradient is used. 

We use a simple model to discuss the possible limitations on the form of the viscous 
stress arising from the requirement that viscous information must travel at a finite speed. 
Contrary to results in the existing literature, the viscous stress tends to be increased rather 
than reduced for the type of fiows we consider in which the angular momentum and angular 
velocity gradients have opposite signs. However, finite propagation effects may result in a 
unique location for the sonic point. 

We found good agreement between the radial fiow and specific angular momentum profiles 
in the inner regions of the one dimensional models and those in the equatorial plane for 



corresponding two dimensional models which may be matched for a range of a between 0.1 
and 10-^. 

1 Introduction 

It has long been realized that the inner regions of accretion disc flows surrounding relatavistic 
objects have different properties from those pertaining to the case when the central object 
can be modelled as a Newtonian point mass. 

In order to discuss such cases the approximation is often made that the flow may be treated 
as Newtonian but with the modification that the gravitational potential may be written in 
standard notation following Paczyhski and Wiita (1980) as 

GM 

where tq is the gravitational radius. In this potential the angular momentum of circular 
orbits has a maximum at Sr^ and so from a particle point of view orbiting material would 
be expected to become unstable and flow towards the centre attaining sonic velocities fairly 
rapidly. Hence a transonic flow is expected in the central regions ( Abramowicz and Zurek, 
1981). 

Models of thick accretion discs with steady transonic flows were constructed by Paczyhski 
(1980) who assumed that matter flowed over the surface of an essentially hydrostatic disc, 
much as a star in a close binary system overflows its Roche lobe. Later work assuming 
the accretion occurred entirely along the equatorial plane was carried out by Paczyhski and 
Abramowicz (1982) and Rozyczka and Muchotrzeb (1982). In this some assumptions are 
usually made so that the problem is effectively reduced to a one dimensionlal one. 

More recently interest has focussed on so called 'slim discs' ( Abramowicz, Czerny, Lasota 
and Szuszkiewicz, 1988). These studies reduce the problem to a one dimensional one, through 
a form of vertical averaging even though the discs are often geometrically thick. Radial advec- 
tion is included and a form of the Shakura and Sunyaev (1973) prescription for the viscosity 
is often used. Modelling of accretion discs on this basis has been carried out more recently by 
Kato, Honma and Matsumoto (1988) and Chen and Taam (1993). It is one of the purposes of 



this paper to construct two dimensional axisymmetric models of the inner region which take 
the vertical structure fully into account in order to compare the properties with those of the 
one dimensional models. We make the simplifying assumptions that the innermost regions 
have constant specific angular momentum and entropy. Apart from reasons of simplicity, the 
justification for this is that when the fiow velocity approaches the sound speed, the infiow or 
advection timescale becomes significantly shorter than the viscous or thermal timescales for 
the region of interest (see Paczyhski and Abramowicz, 1982) so justifying the assumptions. 
We remark that the inner regions of the one dimensional models with viscosity included are 
found to have a very slowly varying specific angular momentum profile. 

Another issue raised recently by Narayan (1992) is the fact that one must contrain ones 
treatment of viscosity so that viscous information does not travel with supersonic velocities 
(see also Chen and Taam, 1993). Narayan (1992) suggests that the viscous stress should 
vanish when the fiow velocity attains the sound speed. This may have important implications 
for models of transonic accretion discs esspecially when sonic velocities are attained outside 

In this paper we develop a simple model of the viscous process which has the propagation 
limitation built in. We find that the viscous stress does not vanish at a sonic point but 
actually tends to be increased if the angular velocity and angular momentum gradients are of 
opposite sign as occurs in the fiows discussed here. However, limitations of finite propagation 
speed may affect issues associated with the uniqueness of the fiow subject to specification of 
appropriate boundary conditions. 

In section 2 we consider steady state accretion tori with constant specific angular momen- 
tum and entropy. We derive an equation for the velocity potential in the special case of zero 
vorticity. In section 3 we discuss the numerical solution of this equation. For the purpose 
of comparison of these models with one dimensional models, we consider such models, which 
include viscosity treated according to the usual Shakura and Sunyaev (1973) a prescription, 
derived from vertically averaged equations in section 4. We argue that specification of the 
mass accretion rate and constant entropy of the fiow coupled with a specification that the 
fiow have Keplerian rotation at some outer boundary radius does not yield a unique location 
for the critical point. 



Further, we discuss, using a simple model, the possible effects of the limitation that viscous 
information should not be transmitted at a speed exceeding the sound velocity. We find that 
if this speed coincides with the sound speed, the lack of uniqueness in determination of the 
critical point discussed above may be removed. 

In section 5 we describe our results for one and two dimensional steady state accretion 
flows. When the intersection point of the sonic surface with the equatorial plane in a two 
dimensional model coincides with the critical point in a one dimensional model and the sound 
velocities match there good agreement is found for the flow parameters. 

Finally in section 6 we discuss our results. 

2 Tori with constant specific angular momentum 

In axisymmetric tori with constant specific angular momentum the velocity u = {ur,u^,Uz) 
in cylindrical coordinates (r, 0, z) is such that / = ru^f) is constant. We assume the polytropic 
equation of state, 

P = /Cp^. (1) 

Here P is the pressure, p is the density, /C and 7 are the polytropic constant and adiabatic 
index, respectively. The sound speed, Cg, is then a simple function of the density and can be 
written in the form 

The equation of motion governing the steady state may be written 

a; X u = -VB, (2) 

where a; is the vorticity and the Bernoulli function 

1 r^ 

2 7-1 

with \1' being the gravitational potential. For this we adopt the 'pseudo Newtonian potential' 
(Paczyhski and Wiita, 1980) given by 

GM 



V^2 _|_ ^2 _ ^^' 
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where G is the gravitational constant, M is a mass of the central object and tq is the 
gravitational radius. Paczyhski and Wiita (1980) emphasized that all relatavistic effects 
which are relevant to the accretion process near the inner edge of a disk can be reproduced 
in the Newtonian formalism assuming the pseudo-Newtonian potential given above. 

From equation (§) it follows that E is constant on stream lines and may be specified to 
vary arbitrarily from stream line to stream line. The fact that it is possible to specify an 
arbitrary function in this way, being associated with a non zero (p component of vorticity, 
results in a degree of arbitrariness in possible steady state solutions. 

2.1 The zero vorticity case 

In this paper we shall consider the case when B is constant. This corresponds to zero vorticity. 
This is the simplest assumption that can be made about B and it is also the one which is most 
likely to result in models that can be approximately described by one dimensional vertically 
averaged ones. 

It must be noted that the ability to specify the Bernoulli function arbitrarily on stream 
lines introduces a high degree of arbitrariness into the two dimensional models. In principle 
one might for example be able to construct models with flow concentrated towards the disc 
surface (Paczyhski, 1980) or towards the equatorial plane (Paczyhski and Abramowicz, 1982) 
but such solutions would have a non zero component of vorticity in the y? direction which 
can only be specified as an entry condition for the flow under the assumption of the lack of 
importance of viscosity we make in the treatment of the inner region. 

When the vorticity is zero, the meridional component of the velocity field may be expressed 
as the gradient of a potential, $, such that 

{ur.Uz) = {d^/dr,d^/dz) . 

For a steady flow, we also have the continuity equation, which for this flow reads 

V(pV<l>) = 0. (3) 

This may be combined with the Bernoulli equation, 

1 c^ /2 

(V$)2 + —^ + ^ + = i3 = constant 

2 7 — 1 2r^ 



to give a single second order partial differential equation for the velocity potential $ : 



v'$ - ^ ("vi^t) + V 



^(V$)^' 



V<l> = 0, 



(4) 



where the total centrifugal plus gravitational potential is 

/2 



^. 



^ + 



2r2' 



In addition the density of the flow can be found from the Bernoulli equation which reads 



P 



:L_l{B-^r-\{V^f 
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(5) 



In order to discuss possible solutions of equation(^, we first consider properties of the 
potential ^t for constant specific angular momentum /. 

The equipotential surfaces have a saddle point or cusp if / is in the range given by Imb < 
I < Ims, where Imb and Ims are the values of the specific angular momentum at the marginally 
bound and marginally stable circular particle orbits respectively ( Abramowicz et al. 1978). 
In order for accretion to take place it is expected that / be in the above range. In order to 
construct a steady state solution we first define a specific angular momentum in the above 
range and use an equipotential surface as a reference surface in order to construct a sonic 
surface. 

2.2 The sonic surface 

The accretion flow is expected to be transonic, starting at small subsonic speed at large radii 
and accelerating to attain the sonic speed along the sonic surface, being the generalization 
of a critical point in one dimension. In general a unique form for the sonic surface cannot 
be constructed in advance of knowing the solution, even if the other bounding surfaces are 
known (see Anderson, 1989). However, construction turns out to be simple in the special 
case when the flow is normal to the surface, (see Anderson, 1989 for reasons why this might 
be a preferred case). A surface of this form may be constructed by first selecting a bounding 
equipotential surface through which the sonic surface passes orthogonally. We suppose this 
to be a zero density surface where the fiow approaches zero velocity, the sound speed being 
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zero. If the equation of the sonic surface with the above properties is r = r{z), on the surface 
we must have 

dr 5$ 9$ 

dz dz dr 

Differentiating this equation with respect to z and using equations (Q) and @, together with 
the fact that on the surface u^ + u1 = c^, we find the following second order differential 
equation for the sonic surface: 

r 



dz^ _ K dr dz dz ) i_ 



l+(f) 2-l(B-*.) 



Here B is constant and equal to the value of \E't on the bounding equipotential surface. 
The sonic surface that passes orthogonally through the chosen bounding zero density surface 
as well as the equatorial plane can be constructed uniquely by solving this equation. The 
zero density surface of the torus, sonic surface, the equator {z = 0) and a bounding surface 
r = constant on which the inflow velocity is small are taken to define the flow domain. We 
assume reflection symmetry in z for the flow below the equator. At the sonic surface we set 
the outflow velocity equal to the sound speed and normal to the surface. 

In this way the parameters defining a solution are the selected constant specific angular 
momentum and the bounding zero density surface . Alternatively one may use the point 
on the equatorial plane through which the sonic surface passes (identified with the critical 
point in one dimensional flows) r = r^,, and the value of GM /{r^c^^) there as the parameters 
defining a model. The solutions can be scaled to give any mass accretion rate because they 
are invariant to an arbitrary scaling of the density. 

3 Numerical method for the two dimensional tori 

In order to solve the equilibrium equation for $ and at the same time find the velocity of the 
fluid, we rewrite equation (^) in the form of a diffusion type equation: 

^ = D(r,^,t)[V(pV<l>)], 



where D{r, z) is a diffusion coefficient, which in general can be a function of r, z and t. The 
form of this can be chosen for numerical convenience subject to the requirement of numerical 
stability. We solved this equation in the elliptic (subsonic) regime of the computational 
grid by solving the diffusion equation as an initial value problem, integrating forward until 
a steady state was achieved. A two stage procedure was used. During the first stage, the 
equations were advanced with p fixed at the value it would have were the flow sonic. Once 
a steady state had been reached, the evolution was continued for the second stage during 
which the density was calculated from equation(^) until a final steady state was attained. It 
is required to specify boundary conditions on the boundary of the computational domain. 
These were that the velocity is sonic and normal to the sonic surface, that the inflow be 
normal at the exterior boundary, that the flow be tangential at the bounding equipotential 
surface and symmetry with respect to reflection in the equatorial plane. The solution is 
undertaken on an equally spaced computational grid with between 40 and 50 grid points in 
the subsonic region of the flow on the equatorial plane. It is straightforward apart from the 
complication that because the density and therefore the sound speed go to zero somewhat 
interior to the bounding surface that would occur if there was hydrostatic equilibrium, there 
must be a supersonic transition slightly interior to this surface. But note that this transition 
occurs at low absolute speeds. The region between this other sonic surface and the bounding 
surface is small and we treated it by not allowing the density to fall below the value taken on 
when the flow is sonic. In this way a transition to a hyperbolic region is prevented. However, 
the location of this new bounding sonic surface could be found in this way and the solution 
continued into the hyperbolic domain by initial value techniques. As this region was small, 
having a negligible effect on the interior, and in reality it would depend on details of the 
proper boundary condition, we did not calculate it in detail. 

Before discussing the solutions for two dimensional tori that were obtained numerically, 
we describe the one dimensional models which correspond to them. In this regard we note 
that the full disk can only be approximated as a constant angular momentum torus in the 
innermost region. Further out there must be material with a higher specific angular momen- 
tum and viscous effects to allow it to accrete. With reference to one dimensional disc models 
we shall argue below that once the two parameters r^, and GM/{r^,c1) are specified and if 



the angular velocity gradient at the critical point corresponds to constant specific angular 
momentum, there is a unique standard " a" viscosity parameter corresponding to the torus. 

4 Basic equations for one dimensional flow 

The basic equations for one dimensional flow with the time dependence retained are the 
equation of motion 

dt dr S Or dr 

Here u is the radial velocity, S is the surface density and P is now the vertically integrated 
pressure which we shall assume to be a function of both S and r. The latter dependence on r 
may come about through the process of vertical averaging (see below) . A local sound speed 
Cs is then defined through 

^ = (^] 



The continuity equation is 



9S (9S _ S djur) 
dt dr r dr 



The specific angular momentum l{r,t) = r^Q{r,t) satisfies 

dl dl 1 dir^Tr^) 

+ U-, — 



dt dr (rS) dr 

Here t^.^ is the rip component of the viscous stress tensor. This is usually taken to be given 
by T.f^ = Tr^o = iy^r{dQ/dr), with u being the kinematic viscosity. For the Shakura and 
Sunyaev (1973) a prescription, we have u = ac^/Qx, ^k being the circular orbit frequency. 

4.1 The question of causality 

In an important paper, Narayan (1992) has pointed out that the conventional formulation 
of viscosity as a diffusion process allows propagation of information at infinite speed and 
so could lead to unphysical results in transonic fiows which attain velocities exceeding the 
maximum velocity of the information carriers which produce the viscosity. Unphysical results 
might arise through unrealizable diffusive communication between different parts of the fluid. 



Using an ilustrative model of the diffusion process with a finite propagation speed, Narayan 
(1992) concludes that the np stress should vanish at a sonic/critical point, where v? = 
Cg, in a one dimensional accretion flow. Here we provide a different but similar example, 
within a Newtonian framework, which uses all the basic equations. The full set of equations 
is of hyperbolic type and so information propagates at finite speeds. The example shows 
that for flows in which the angular velocity and specific angular momentum gradients are 
of opposite sign, the viscous stress is, if anything increased, rather than reduced at the 
critical point. Furthermore, when the specific angular momentum gradient is very small 
when the radial velocities are significant, such as might be expected for the flows considered 
here, the corrections due to finite propagation effects become small so that the conventional 
Shakura and Sunyaev (1973) treatment should be reasonable. However, constraints arising 
from causality considerations may be important when the question of the uniqueness of steady 
state solutions is discussed. When the angular velocity gradient is very large, the conclusions 
derived from this simplified model are similar to those of Narayan (1992). 

The basic equations of our simplified model are equations (|), (^ and (H) together with 
an equation for Tn^ of the form 

OTrip OTj-in vripO '^np) /ri\ 

^ + "^ = r • ^'^ 

This states that the stress relaxes locally to its equilibrium value on a relaxation timescale r 

which may be an arbitrary function of the state variables but not their gradients. We stress 

that this formalism has been introduced to ensure that information propagates at a finite 

speed in the system rather than to provide a realistic formulation of viscosity. It is thus 

adequate to demonstrate that zero viscous stress at the sonic point cannot be inferred from 

any limitation in the propagation speed of information. Equations (^, (^, (||) and (|) form 

a system of four simultaneous first order partial differential equations which may be written 

in the form: 

f^ + E^.-^ + 5^ = 0.(^ = l,2,3,4) (10) 

i=i 

Here t/^ = (S,M,/,r,^), Si = Jlu/r, 82 = ^ (f )^ - rfi^ + ^, 5-3 = -2r,^/S, 5-4 = 
{Tr^ + 2z/S/r^^)r^^. The diagonal elements of the matrix (A) are all equal to u, Ayi = 
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S, A21 = c^/S, A34 = —r/T,, A43 = —uT./^rr), with all other elements being zero. The 
system is hyperbolic if, as for our set of equations, the eigenvalues of A are all real. The 
characteristics are rays with propagation speeds given by these eigenvalues. In our case these 
are the two sonic speeds u + Cs,u — Cs and the two viscous speeds u + c^, u — c^, with 



c„ = v'^/t, being the speed of propagation of viscous information. Therefore information 
cannot propagate faster than a finite speed in the model system. But note that the viscous 
speed, which we assume to be comparable to the sound speed, becomes infinite when the 
relaxation time is zero. When c^ = Cg., we have a = QkT- 

4.2 The steady state 

We now examine the steady state in which ^ = 0. The continuity equation gives on integra- 
tion 

2nEru = -ci, (11) 

Ci being the constant accretion rate. Equation (j^) gives on integration 

ci/(r) + 2'Kr'^Tr^ = C1C2, (12) 

where C2 is the constant rate of flow of angular momentum through a circle of radius r. 
Equation (||) gives 

Trip = uT^r{dVL/dr) — tu—^. (13) 



From the above two equations it follows, after differentiating the first and using the result to 

eliminate -^, that: 

u^r{{dnidr)-{u'/{cy)){dl/dr)) 
^''^~ (l-(2izr/r)) ■ ^ ^ 

This is the Shakura and Sunyaev (1973) viscosity prescription modified by additional terms 
involving the radial velocity. We first remark that for u ~ — Cs, and r < fi^^, \uTJr\ < H/r, H 



being the disk thickness. Thus for most purposes the denominator in equation (|1^) may be 
approximated as unity (it could also be absorbed by a redefinition of a incorporating radial 
dependence) . 

We note that corrections arising from propagation effects are small if \dl/dr\ is small such 
the flow has almost constant specufic angular momentum even when u"^ = c^. Furthermore 
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the viscous stress is increased if the gradients of angular velocity and angular momentum 

are of opposite sign. But note that an interesting situation arises at a sonic point in the 

physically reasonable case when c^ = Cs. Equation ([I^ ) may be written 

z/Er ((1 - (MV( cg))(dn/dr)- {2nu'') / {rcD) 
;i - {2uT/r)) 

From this it follows that at a sonic point: 



^ ^ '^^' u-^ ~ I"- / I'-sjn"-"/"-' ; ~ i-^""- ;/ y ^-s)) q^n 



2z/Sfi , . 

-•■- = -Jl^mry (16) 

We see that the stress is independent of the angular velocity gradient and takes on the same 

form as it would do in the case of constant specific angular momentum. Equation(|IBp shows 

that the angular momentum flux is then fixed at the sonic point by the magnitude of Vt. 

This is in contrast to the situation when the viscosity prescription with m = is used. 
Then the angular momentum flux may be adjusted by altering the angular velocity gradient 
at the sonic point, albeit by a small amount if a is small. Within this prescription, when 
dQ/dr = —2Q/r, the angular momentum flux is also fixed by Q at the sonic point and the 
standard a prescription is essentially unmodified by finite propagation effects. We shall refer 
to this condition at the sonic point as the natural boundary condition. 

In fact equation (|15|) implies the existence of a critical point at the sonic point and that 
equation (|r^) should be satisfied there to avoid a singularity and that indeed freedom to 
specify the angular velocity gradient should be lost. However, the technicalities of additional 
critical point analysis ( scarcely justified by the sophistication of the model ) can be avoided 
by noting that for our solutions, when the natural boundary condition is used, equation (^) is 
essentially satisfied in any case and the angular momentum profiles are so fiat that limitations 
due to finite propagation effects may be ignored. Note too that this is reasonable on physical 
grounds because viscous effects are small in any case when the angular momentum profile is 
fiat. 

For the simple model discussed here, we have seen that when the viscous propagation 
speed is equal to the sound speed, freedom to adjust the angular momentum flux by varying 
the angular velocity gradient at the critical point is lost. This flux is always the same as when 
the natural boundary condition is used. We shall see that this is important for the issue of 
uniqueness of steady state solutions. 
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When the mass flow rate, a and the enropy of an adiabatic steady state solution are 
specified and the natural boundary condition is used, the solution is apparantly unique. 
However, if there is freedom to adjust the angular momentum flux by varying the angular 
velocity gradient at the critical point, as might arise if c^ > Cg, the solution is no longer 
unique (see below). 

Finally, we comment that when the angular velocity gradient is very steep, we have 
approximately 



_ z/Sr (1 - ^) jdQ/dr) 

This is of a similar form to that given by Narayan (1992). 

In the numerical work presented below, we adopt the standard Shakura and Sunyaev 
(1973) form: r^p = Tr^po, noting that we expect modifications due to finite propagation effects 
to be negligible, when the natural boundary condition is used. 

For the steady state, the equation of motion reads: 

du 1 dP dm f 

The vertically integrated pressure may be taken to obey a polytropic equation of state P = 
^^i+i/n^ n being the polytropic index. The polytropic factor K may be taken to be constant, 
or a function of r to reflect vertical averaging. This dependence can be found from the scalings 
for an adiabatic equation of state: P oc Hp'^, S oc pH, and Vt^H"^ oc c^ oc p^~^, with VLk = 



jGM/{r{r — TgY) being the Keplerian angular velocity for circular orbits. Eliminating H 
and p from these gives 

2(7-1) 

Kq being constant and 1 + 1/n = (87 — l)/(7 + 1). Combining the equation of motion with 

the continuity equation and the above equation of state gives 

I cg\ du _ GM cl ( nr dK\ 2 

y u ) dr (r — tgY r \ (n + 1)K dr ) 

This exhibits a critical point when the flow speed reaches the sound speed. At such a point 

the right hand side of equation ([T9|) must vanish. The nature of the critical point when 

K is constant has been discussed by Papaloizou and Szuszkiewicz (1993) who find that for 

physically acceptable solutions it is of saddle type. 
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4.3 Dimensionless form of the equations 

In order to discuss solutions of the steady state equations further, we adopt the usual "a" 
viscosity: 

Tr^ = Tr^o = uJ]r{dQ/dr), 

with u = ac^/Qx- We use this together with equations (pTf ), ( [T2|) and ([19|) . 

We find it convenient to adopt as units of radius and velocity r* and m* being the radius 
at the critical point and the absolute magnitude of the velocity velocity there, respectively. 
Then we define the dimensionless variables v = —uju^, x = r/r^,, u = Qr^/u^,, ujk = ^Kf^^ju^, 
and xg = tg/t*- Note that because we are dealing with inflow v is positive. The square of the 
dimensionless sound speed is c^ = i^/(fa;)^'", with K = K{x)/K(l). In addition we define 
C2 = C2/(u*r*), and m = GM/{r\u1) as the Bondi parameter. 

It is possible to derive a single nonlinear second order differential equation for u from the 
above equations which may be written 



^=lLi-Vn- ^ \^^^n ^ f 1 _ ^^Vo^^a; = (20) 

ni\ rf.i/nyi+2/n J ^3, (^x-xgY fV^aj^^V mK dx ) ' 



where, ni = 1 + 1/n and 



yn. ^ '^ ^_ (21) 

ujkX^/^~'^{c2 — ojx"^) dx 



At the critical point where x = v = K = 1, equation ( pOD requires that the dimensionless 
angular velocity oj be given by 



uo* 



m n I dK 

1 + 



\ (1 - xgY (n + l) \ dx ) ^^^ 



Equation (|2lD then gives 

_ ( a duj\ , , 

C2= — — +UJ,. (22) 

\UKdx)^^^ 

We have also found it convenient to work in terms of a quantity cn related to 62 through 

C2 = jT^ + ^*- (23) 

Thus when the natural boundary condition is used ctv = 1- 
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4.4 Parameters specifying a solution 

Equation (^) after use of (0) constitutes a second order differential equation for uj. After 
specification of the Bondi parameter m,XG, ol and C2, both lo and its first derivative are 
specified at the critical point. For the critical points considered here these conditions specify 
a unique solution that is subsonic for x > 1 ( Papaloizou and Szuszkiewicz, 1993). In order 
to be physically acceptable this solution must match onto a Keplerian or near Keplerian 
disk at some specified large radius. In general this requires that one of the parameters such 
as a be adjusted as an eigenvalue in order to fulfill this requirement. In this way we see 
that specification of m, xq and C2 determines a value of a required to satisfy the boundary 
conditions. 

However the parameter xq which can be regarded as specifying the location of the critical 
point may also be varied, other parameters being kept fixed, so determining this quantity as 
a function of a. 

On physical grounds we expect that for adiabatic solutions of the type we are considering, 
we should be able to specify the accretion rate and the entropy with some degree of arbitrari- 
ness. The solutions can always be adjusted to give any accretion rate by scaling the constant 
Ci, while the Bondi parameter m may be used to adjust the entropy. If these are determined 
in this way, the above discussion shows that when 62 is fixed the location of the critical point 
is then determined as a function of a. 

On dynamical grounds we expect the solutions considered here to be of nearly constant 
specific angular momentum when there are significant velocities and the discussion of the 
requirement of a finite propagation speed of viscous information given above then suggests 
that 62 be chosen to satisfy the natural boundary condition such that 

Under these conditions C2 is no longer free and we expect that for a given accretion rate, 
entropy and a, the solution together with the location of the critical point are specified 
uniquely at least locally in parameter space. However, if we allow freedom in the specification 
of C2 this is no longer precisely so. For a given m and xg there is a possible range of a 
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corresponding to the range in 62 or equivalently for a given a there is a possible variation in 
Xq or the location of the critical point. 

Variation of 62 causes the angular velocity gradient to differ from that appropriate to 
constant specific angular momentum at the critical point. Because the dynamics requires 
that the angular momentum profile does not differ much from a constant, this situation is 
restored in a fairly narrow boundary layer when the boundary condition for 62 differs from 
the natural one. In practice, for a fixed m, the variation in a allowed for a given xq-, or 
the variation in xq found for a given a consequent on varying 62 by a maximum reasonable 
amount, is found to be small. 

Adopting the natural boundary condition, the above discussion indicates that specification 
of m and a results in a unique critical point location given through xq. 

Such a solution can then be matched to a two dimensional solution with constant specific 
angular momentum which has the same value of m and which has a sonic surface which passes 
through the equatorial plane at the same location as the critical point in the one dimensional 
case. 

4.5 Numerical methods for the one dimensional case 

Numerical solutions in the one dimensional case have been generated by two methods. In 
the first, the domain between the critical point and the outer boundary was divided into A^ 
equally spaced zones, with values of A^ up to 18000 having been used. Equations ( ^0|) and (21) 



were then written in a centered finite difference form. On specification of CAr,m,XG and a, 
the solution can be found by stepping outwards from the critical point. However, as discussed 
above, the condition that the angular velocity match on to the Keplerian value at some outer 
boundary radius, in general cannot be satisfied unless a takes on a specific characteristic 
value. This method was most satisfactory when the critical point was inside Sr^ and vti was 
not too large. When m is large, the disk thickness is very much less than the radius, and 
the equations are stiff. This has the effect that with the above method it is awkward to 
extend the solution to very large outer boundary radii. However, the characteristic value of 
a corresponding to given xq and ca? is very accurately determined. 
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An alternative approach is to use a relaxation method (Papaloizou and Szuszkiewicz 1993) 
such as for example solving the diffusion type equation derived from equation (|20|) 



where P is a diffusion coefficient which can be a function of x and t and which can be chosen 
for numerical convenience. 

This can be solved in a computational domain extending from the critical point to the 
outer boundary. In this case uj is specified at the domain end points. For fixed a, m and 
xg this is enough to specify the solution but then 62 is determined in the course of solution 
through equation (^11) . This has the disadvantage that in order for a solution to exist, once 
Xg and m have been specified, a must be in the range given through the maximum permitted 
variation of 62 (see below). This method which does not have the stiffness difficulties of the 
first was found to work best when the critical point was outside Src. Some solutions obtained 
by the first method were checked using this method. 

5 Numerical Results 

We have constructed two dimensional model tori with different values of the constant specific 
angular momentum and vertical thickness. Each torus is characterized by a value of m and 
Xg given in Table 1. For each of these models the outer boundary was taken at r = Atg 
which was far enough out for the inflow velocity to be small there, making it seem reasonable 
that the disk can be regarded as essentially static beyond. 

The results we describe below are qualititatively similar for all of the models we examined 
so we shall focus on two examples for illustrative purposes. We show the velocity and density 
structure appropriate to models 3 and 12 in figures 1 and 2 respectively. Model 3 represents a 
moderately thick disc. The inner edge of this disc, where the radial inflow reaches the sound 
speed, is located at the radius ~ 2.5rG'. The maximum density occurs around 3.5rG'. Model 
12 is an example of a genuinely thick accretion disc. Its inner edge is close to 2.05rG and the 
density maximum is outside Atg- 
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For each of the two dimensional models we have constructed a one dimensional viscous 
equivalent with the same values of m and xq- As argued above if cn, which controls the 
angular velocity gradient at the critical point is fixed, there is a unique value of a for which 
the angular velocity takes on the Keplerian value at some outer radius, here taken to be at 
IOtg- We denote the value of a corresponding to cat = 1, the natural boundary condition, by 
amax- In addition we denote the value of a corresponding to cat = or zero angular velocity 
gradient at the critical point by amin- 

The range between amax and amin is spanned continuously as c^ is varied from 1 to 0. 
Values of c^ less than correspond to the situation where the angular velocity gradient turns 
over. If variation in c^ is permitted this results in a range of possible xq for fixed m and a, 
(see the above discussion of uniqueness). We remark that amax and amin do not normally 
differ by more than twenty percent. In addition, other parameters being fixed, these values 
are somewhat sensitive to the location of the outer boundary and tend to decrease as this 
moves outwards. The sensitivity is greater for smaller values of m. The angular momentum 
profile interior to Atq is very fiat with / oc r~'^'°'^. This is fully consistent with the idea that 
viscous effects are weak in this region and justifies the approximation of constant specific 
angular momentum. The one dimensional models do not depend very much on whether 
vertical averaging is taken into account in the equation of state {K = K{x)) or a straight 
forward polytropic equation of state {K = constant) is used. 

Values of amax and amin for some values of xq and m are given for the former case in 
table 2 and the latter case in table 3. 

It is seen that a range of 7.6 x 10^® < « < 1.7 x 10^^ can be spanned by varying m and xq- 
The specific angular momentum for one dimensional models ( with vertical averaging taken 
into account in the equation of state ) and two dimensional models are compared in figures 
3 and 4. In figure 3 we show the distributions of specific angular momentum in the model 3 
torus and in the corresponding one dimensional vertically integrated models with cat equal 
to zero and one respectively. For cjv = a very narrow, but clearly pronunced boundary 
layer exists. In this layer, the specific angular momentum rises steeply, eventually attaining 
an almost constant value. In figure 4 we plot the same curves but for model 12. 

The radial velocities in the one dimensional models are not very sensitive to whether 



Table 1: The various two dimensional model calculations: m is the Bondi parameter measured 
where the sonic surface intersects the equatorial plane at r = r*. The equipotential surfaces 
have a cusp at r = rcTeq, h is the ratio of the vertical distance between the equatorial plane 
and the zero density surface of the torus, to the radius and xq = tg/t^:. In addition we 
show Umax and amm? the values of a appropriate to one dimensional models, with vertical 
averaging taken into account in the equation of state, with the same xq and m. The former 
value corresponds to the natural boundary condition while the latter corresponds to zero 
angular velocity gradient at the critical point. For each one dimensional model, the exterior 
boundary was taken at lOrc. 

Model m Xq r^q h a^ax Oimin 

1 



"" -^Cr ' eq '" '-^max '^mm 

302 0.34 2.9 0.1 6.03 x lO^^ 2.98 x lO^^ 

83 0.39 2.9 0.2 1.26x10-2 1.07 x 10"^ 

3 75 0.39 2.7 0.2 1.40 x 10"^ 1.17x10^2 

4 21 0.42 2.7 0.4 1.81 x 10"^ 1.50 x 10"^ 

5 63 0.43 2.4 0.2 3.65 x lO"''^ 3.48 x lO^^ 

6 16 0.43 2.7 0.5 1.58 x lO^^ 1.34 x 10"^ 

7 18 0.45 2.4 0.4 9.01 x lO^^ 8.18 x lO^^ 

8 13 0.45 2.4 0.5 9.94 x lO^^ 8.97 x lO"''^ 
55 0.46 2.2 0.2 5.67 x 10"^ 5.63 x 10^^ 
51 0.48 2.1 0.2 2.59 x 10"^ 2.58 x 10"^ 

11 15 0.48 2.2 0.4 3.50x10-''^ 3.38x10^3 

12 11 0.49 2.2 0.5 4.02 x lO^^ 3.86 x lO^^ 



8 
9 

10 
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Table 2: Values of Umax and amin-, for a sample of one dimensional models with various xq 
and m for the case when vertical averaging was taken into account in the equation of state. 
In each case 7 = 4/3 and the external radius was taken at lOrc. 



m Xq O^max ^min 

10 0.32 1.67 X 10^1 6.24 x 10^^ 

100 0.35 6.51 X 10-2 3.39 x lO'^ 

10 0.35 1.00 X 10-1 5.08 x lO^^ 

100 0.40 7.45 X 10-3 6.71 x lO'^ 

10 0.40 4.29 X 10-2 3.00 x 10"^ 

100 0.45 2.71 X 10-^ 2.71 x 10"^ 

10 0.45 1.44 X 10-2 1.25 x 10-^ 

100 0.50 7.61 X 10-6 7.61 x 10-<^ 

10 0.50 3.18 X 10-3 3.08 x 10-^ 
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Table 3: Values of amax and amin, for a sample of one dimensional models with various xg 
and m for the case when vertical averaging was not taken into account in the equation of 
state. In each case 7 = 4/3 and the external radius was taken at lOrc- 

m Xq O^max ^min 

10 0.32 1.23 X 10^^ 5.51 x 10^^ 

100 0.35 5.96 X 10-^ 3.26 x lO'^ 

10 0.35 7.88 X 10-2 4.46 x 10^^ 

100 0.40 8.01 X 10-3 7.23 x lO'^ 

10 0.40 3.55 X 10-2 2.64 x 10-^ 

100 0.45 4.71 X 10-^ 4.70 x 10"^ 

10 0.45 1.31x10-2 1.17x10-2 

100 0.50 2.41 X 10-5 2.33 x 10-^ 

10 0.50 3.54 X 10-3 3.43 x 10-^ 
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Cat = or Cat = 1 is used. This is illustrated in figures 5 and 6 in which the radial velocities 
are plotted for the one dimensional models corresponding to models 3 and 12 respectively for 
a small range of radii near to the critical point. 

In figures 7 and 8 radial velocity profiles in one and two dimensional cases are compared 
for models 3 and 12 respectively. For the two dimensional tori, we plot the radial velocity 
measured on the equatorial plane. We plot the radial velocities for the corresponding one 
dimensional models with and without vertical averaging being taken into account in the 
equation of state. There is good agreement between all three of the radial velocity profiles 
for each model. This indicates that the results are not very sensitive to the precise mode of 
vertical averaging. 

6 Discussion 

In this paper we constructed steady state non static accretion tori with constant specific 
angular momentum and entropy. In general there is some degree of arbitrariness in this 
procedure because the Bernoulli integral may be chosen to vary in an arbitrary way on 
stream lines. This corresponds to differring specifications for the if component of vorticity 
which needs to be specified as an entry condition on the fiow. For simplicity we have restricted 
attention to the case of zero vorticity. In this case we derived an equation for the velocity 
potential and solved it to give solutions appropriate to the inner transonic region of an 
accretion disc. 

For the purpose of comparison of these models with models obtained in a one dimensional 
approximation, we considered such models, including viscosity treated according to the usual 
Shakura and Sunyaev (1973) a prescription. 

An important issue is whether specification of the mass accretion rate, entropy ( or sound 
speed at the critical point) and the constraint that the material have the specific angular 
momentum appropriate to a Keplerian disc at some specified outer boundary radius is enough 
to fix the location of the critical point for a given value of a. 

We found that if the usual a viscosity prescription is used in which the viscous stress 
is proportional to the angular velocity gradient, the location of the critical point is not 
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determined uniquely, although the possible spread in locations might be small. Note too 
that previous work by Abramowicz et al (1988) used a viscosity prescription in which the 
viscous stress was simply proportional to the pressure. In this case one expects a unique 
determination for the location of the critical point (at least locally in parameter space.) The 
extra freedom in our case arises from the possibility of varying the viscous stress at the critical 
point through changing the angular velocity gradient. 

In order to investigate this situation further, we discussed, using a simple model, the 
possible effects of the limitation that viscous information should not be transmitted at a speed 
exceeding the sound velocity (Narayan, 1992). We found that contrary to some previous work, 
the viscous stress does not vanish at the critical point but that it tends to be increased if as 
in our models, the angular momentum and angular velocity gradients have opposite sign. In 
addition models such as the ones we constructed with very fiat specific angular momentum 
profiles throughout ( ie with no boundary layer) are barely affected by this propagation 
constraint. 

However, if the speed associated with the propagation of viscous information coincides 
with the sound speed, the ability to change the viscous stress through varying the angular 
velocity gradient at the critical point together with the lack of uniqueness in determining its 
location may be removed ( at least in our simple model). 

In general we found good agreement between fiow parameters determined from the one and 
two dimensional models. When the intersection point of the sonic surface with the equatorial 
plane in a two dimensional model coincides with the critical point in a one dimensional model 
and the sound velocities match there good agreement is found for the radial velocity profile 
measured on the equatorial plane in the two dimensional case and the radial velocity profile 
found in the one dimensional case. 
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Figure captions 

Figure 1: The density structure (gray scale levels) and the velocity field (arrows) for the 
model 3 torus with m = 75 and xq = 0.39. The dashed curve is the sonic surface. 
Figuere 2: As in figure 1 but for model 12 which has with m = 11 and xq = 0.49. 



Figure 3: The specific angular momentum distribution, given in units of \fGMrG^ for 
the model 3 torus (dashed line) and corresponding one dimensional models (dot-dashed and 
solid lines). The dot-dashed line represents the solution with the natural boundary condition 
ctv = 1 and the solid line the solution with cat = 0. 

Figure 4: As in figure 3 but for model 12. 

Figure 5: The magnified region near the critical point showing the comparison between 
the radial velocities in the one dimensional models obtained with cat = 1 (dot-dashed curve) 
and Cjv = (solid curve) for model 3. 

Figure 6: As in figure 5 but for model 12. 

Figure 7: The radial velocity in units of the sound speed at the critical point is plotted 
for the model 3 torus (dashed curve), for the one dimensional model with vertical averaging 
taken into account in the equation of state (solid curve) and for the two dimensional polytrope 
(dot-dashed curve): 

Figure 8: As in figure 7 but for model 12. 
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